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ICASE Fluid Mechanics 


Due to increasing research being conducted at ICASE in the field of fluid mechanics, 
future ICASE reports in this area of research will be printed with a green cover. Applied 
and numerical mathematics reports will have the familiar blue cover, while computer science 
reports will have yellow covers. In all other aspects the reports will remain the same; in 
particular, they will continue to be submitted to the appropriate journals or conferences for 
formal publication. 
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ANALYTICAL DESCRIPTION OF THE BREAKUP 
OF LIQUID JETS IN AIR 


Demetrios T. Papageorgiov} 
Department of Mathematics and 
Center for Applied Mathematics and Statistics 
New Jersey Institute of Technology 
Newark, New Jersey 07102 


ABSTRACT 

A viscous or inviscid cylindrical jet with surface tension in a vacuum tends to pinch due 
to the mechanism of capillary instability. We construct similarity solutions which describe 
this phenomenon as a critical time is encountered, for two physically distinct cases: (i) 
Inviscid jets governed by the Euler equations, (ii) highly viscous jets governed by the Stokes 
equations. In both cases the only assumption imposed is that at the time of pinching the 
jet shape has a radial length scale which is smaller than the axial length scale. For the 
inviscid case, we show that our solution corresponds exactly to one member of the one- 
parameter family of solutions obtained from slender jet theories and the shape of the jet is 
locally concave at breakup. For highly viscous jets our theory predicts local shapes which 
are monotonic increasing or decreasing indicating the formation of a mother drop connected 
to the jet by a thin fluid tube. This qualitative behavior is in complete agreement with both 
direct numerical simulations and experimental observations. 


lr This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-19480 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681. 
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1 Introduction 


It is well known (Rayleigh (1878)), that a circular jet of finite radius with a surface which 
supports surface tension is linearly unstable to a long- wave capillary instability. Any small 
perturbations with wavelengths larger than the jet radius grow exponentially. Furthermore, 
linear theory predicts a maximally growing wave and hence a dominant length scale for 
the instability. According to linear theory, then, the dominant wavelength is approximately 
9a where a is the unperturbed jet radius (this result is independent of the surface tension 
coefficient). Experiments indicate that the instability can lead to breakup, or pinching, of the 
jet into drops. Clearly the pinching phenomenon is nonlinear since the initial disturbance has 
to grow to amplitudes of the order of the unperturbed jet radius. Linear theory, however, 
does well in the qualitative prediction of breakup times, for instance, by employment of 
empirical arguments such as e-fold amplification of perturbations. In many applications the 
shape and jet velocities at breakup are useful but cannot be obtained from linear theory. 
This classical problem has been studied extensively; experiments have been carried out by 
Donnelly and Glaberson (1966), Goedde and Yuen (1970) and more recently Chaudhary and 
Maxworthy (1980a,b). Weakly nonlinear theories (see below) have been carried out by Yuen 
(1968) and later by Chaudhary and Redekopp (1980). A review of the subject can be found 
in Bogy (1979) while recent simulations using boundary integral techniques are described in 
Mansour and Lundgren (1990). 

Viscosity dominated flows form a separate but parallel field. Tomotika (1936) consid- 
ered the linear stability of a stationary cylindrical thread of viscous fluid surrounded by a 
second viscous fluid with surface tension acting at the interface. Qualitatively, the stability 
results are similar to inviscid studies with a maximally growing wave of the order of the 
unperturbed thread radius. A more complete theory, including the effects of non-uniform jet 
velocities can be found in Chandrasekhar (1961), where it is shown that capillary instability 
provides linearly growing waves that scale on the jet radius. Recently, Tjahjadi, Stone and 
Ottino (1992), have undertaken an experimental and numerical study of breakup of viscous 
cylindrical threads of one fluid in another. The experiments and the computations show that 
at the time of breakup the jet tends to form larger mother drops joined to smaller satellite 
drops by thin slender tubes. It is precisely this regime that we can describe theoretically 
with excellent qualitative agreement with both computations and experiments. Besides their 
intrinsic interest, such solutions can be useful in providing initial conditions for continuation 
of numerical solutions just before the change in topology necessitated by the pinching. The 
analytical description of local structures is also useful in the determination of the effect of 
additional physicochemical effects such as surface active agents or electrical forces on the 
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pinching process. 

The present approach is a fully nonlinear one with interfacial deflections as large as the 
undisturbed jet radius, as opposed to most previous weakly nonlinear studies. Chaudhary 
and Redekopp (1980) (see also Yuen (1968)) consider two low amplitude (asymptotically 
small but not infinitesimally so) initial perturbations, a fundamental and a harmonic. These 
are followed up to cubic order in the initial small amplitude, which is the first stage when 
the modes interact nonlinearly to produce an amplitude equation. The methodology is that 
of the Stuart (1960), Watson (1960) classical weakly nonlinear theory, even though the jet 
problem always has a band of unstable waves which does not become monochromatic as a 
flow parameter (e.g. capillary number here) is varied. The results are therefore expected to 
be valid for sufficiently small times, but as Chaudhary and Redekopp indicate qualitative 
features of the experiments are reproduced at times which seem to be beyond the validity 
of the theory. Fully nonlinear theories allow the interfacial amplitude to be as large as the 
unperturbed jet radius, a situation which is essential in the description of breakup. This 
usually means that the problem should be addressed numerically, and it is the objective of 
this work to present an analytical alternative to breakup. A fully nonlinear theory of jet 
breakup in a vacuum was developed by Ting and Keller (1991), referred to as TK. This 
is done by use of slender jet theory (see below) which was extended by Papageorgiou and 
Orellana (1993), referred to as PO, to describe breakup of jets of one fluid into another 
with different density, with or without surface tension and inclusion of Kelvin-Helmholtz 
instability. In both cases, a simplified set of evolution equations involving time and the axial 
coordinate alone, are derived by an asymptotic expansion procedure which uses the ratio 
between undisturbed jet radius to characteristic axial length scale as a small parameter. As 
noted by Papageorgiou and Orellana, such an ansatz can be applied to flows which have 
initial conditions characterized by a long axial length scale, for instance. This implies that 
the slender jet theory preserves such geometry throughout the evolution and at the time of 
pinching the distance between successive pinch points is large compared to the undisturbed 
jet radius. Qualitatively, therefore, the slender jet theory predicts drops with radii larger 
than the undisturbed jet radius. The strength of the slender jet theory, however, is in its 
prediction of local structures at the time of pinching. As found in TK and extended to other 
flows by PO, a one-parameter family of similarity solutions is possible. Different members 
of this family correspond to local geometrical shapes of the jet which are concave, conical 
and convex respectively (see PO for representative numerical solutions). In this work we 
show that if the slender jet theory is not adopted throughout the evolution, the jet can 
pinch according to a single member of the one-parameter slender jet family, corresponding 
to concave shapes. Solutions with zero surface tension also follow by simply discarding one 
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term from the similarity equations. As shown in PO, the slender jet equations with zero 
surface tension constitute a system of two nonlinear hyperbolic conservation laws with the 
final outcome of the initial value problem being the formation of a shock after a finite time. 
Clearly the development of short axial length scales is inconsistent with the fundamental 
assumptions of the theory. The direct theory presented here, however, does not assume 
the slender jet approximation throughout the evolution and so does not suffer from the 
restrictions outlined above. Consequently, a pinching solution for zero surface tension can 
be given for a flow which is stable according to linear theory. The shapes at breakup are 
monotonic indicating the possibility of formation of mother drops joined to satellite drops by 
thin fluid tubes. Such shapes are seen in experimental studies of inviscid jets (see Donnelly 
and Glaberson (1966)). 

We emphasize that slender jet theories are powerful methods in the description of local 
structures in pinching phenomena. In capillary instability phenomena in three dimensions, 
slender jet theories provide simplified sets of evolution equations which can be analyzed for 
breakup. Since at breakup of cylindrical jets the slender jet assumptions are valid it is not 
surprising to find that a description of breakup emerges. Physically, the disturbance must be 
characterized by large wavelengths initially. The present approach corresponds to a different 
physical set-up as we describe next. If the jet is allowed to undergo natural growth without 
any external forcings and linear waves are at first important, then linear theory can be used 
to predict a maximally growing wave. As mentioned above, this wave has wavelength of the 
order of the undisturbed jet radius and so when the evolution enters the nonlinear regime, 
the interfacial waves have axial length scales comparable to the undisturbed jet radius. 
The slender jet approximation may not be appropriate, therefore, for the total duration of 
the evolution to breakup, even though it is appropriate locally at breakup. The analysis 
presented here, assumes that near the singular time the jet radius is asymptotically smaller 
than the axial length scale, a condition which may be inappropriate for the past history of 
the evolution. 

The article is organized as follows. Section 2 deals with the breakup of inviscid jets 
governed by the axially symmetric Euler equations. A system of differential equations is de- 
rived that describes the local structure of the shape and the potential of the jet at breakup. 
Numerical solutions are given and a correspondence is made with the slender jet results of 
PO, both in the presence and absence of surface tension. In Section 3 we analyze viscous 
axisymmetric jets. The equations studied are the Stokes equations (creeping motion corre- 
sponding to small Reynolds numbers) and a differential equation describing the jet shape at 
breakup emerges. The analysis leading to this equation is more involved than the inviscid 
case due to the higher order of the field equation and the additional stress boundary condi- 
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tions. Numerical and analytical aspects of the equation are also given. Finally Section 4 is 
devoted to discussion and comparison with experiments and simulations. 


2 The theory for inviscid jets 


The equations governing the fluid motion can be written in terms of a velocity poten- 
tial <f>(t,r,z) which is independent of the azimuthal angle 9 for axisymmetric flows. A 
system of cylindrical polar coordinates (r, 9, z) is used with corresponding velocity vector 
u = (it, v, w') = ^ (f). Using this notation, </» satisfies Laplace’s equation while at the free 
surface two conditions are specified, a kinematic condition and a normal stresses balance 
which gives the pressure jump across the interface. When the latter condition is used in the 
2 -momentum equation the usual Bernoulli equation arises. Without loss of generality we 
give the system for a static undisturbed jet, noting that any background axial velocity can 
be removed by a Galilean transformation. The equations with boundary conditions are 

<t>rr + ~<f>T + 4>zz = 0, (2-la) 


On r = S(t, z), 


<j>r — &t + <!>zSzi 

* + !«? + <!) — g -A. + ^)(> + «)-'’■ 


( 2 . 16 ) 

(2.1c) 


A final condition which completely specifies the problem is regularity of <j) at r — 0. Equations 
(2.1a-c) are in non-dimensional form with the non-dimensional group containing the surface 
tension coefficient (the Weber number) scaled to unity, in anticipation of this canonical 
rescaling in the final similarity equations. 

In general (2.1a-c) must be addressed numerically. We construct pinching solutions of 
(2.1a-c) by assuming the following ansatz as a singular time t s is approached from below: 


r = (t,-t)°y , z = , S(t,z) = {t a -t) a f(0, (2.2a) 


(f>{t,r,z) = (t, - t)' r x(Uj/)0- 


According to these transformations time and space derivatives become 

8_ d_ d 

dt dt (t s — t)d£' 



(ts - ty a 


8y 


d 

dz 


(ts - tr* 


d 

df 


(2.2 b) 


(2.2 c) 
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Physically, we expect the jet radius to be much smaller than the axial length scale at breakup. 
This implies that r « z which in turn yields the condition 

a> (3. (2.3) 

Inequality (2.3) is the only necessary assumption for the theory to be consistent. 

Substitution of the ansatz (2.2a-c) into the Laplace equation (2.1a) indicates that since 
a > (3, the dominant terms come from the ?’-derivatives. The leading order contribution 
which is regular at the jet axis y = 0 is independent of y , and we must go to a higher order 
to obtain a non-trivial balance. This suggests seeking a solution by expanding \ in (2.2b) 
in appropriate powers of ( t s — t) as follows: 

<t> = (t s - ty (xo(y,£) + (ts - t) 0 xi(y>O + (ts - t) a x2(y,0 + •••)• (2.4) 

Substitution of (2.4) into the Laplace equation (2.1a) gives 

(t s — ty 2a + ~X°yj + (ts — ty + P 2a |xis/3/ + ~Xly^ 

(t 3 - ty~ 20 x o« + . • • = 0. (2.5) 

It can be seen from (2.5) that the first stage where a regular ^-dependent solution is possible, 
arises if the Xi terms are in balance with the xott term. This yields 

a = 2/3. (2.6) 

The solutions for xo and Xi which are regular at y = 0 are 

Xo = Xi = — + B(0i (2.7 o, b) 

where A and B are unknown functions of £. 

Next we consider the kinematic condition (2.1b). Substitution of the ansatz (2.2) and 
(2.4) along with the solutions (2.7a,b) yields the equation 

((, - (- 2 / a «) = - a (t. - tr~'f 

+ ((,-()’ + “- w -4 ( /j + ... . (2.8) 

A few comments about (2.8) are in order. First, the two terms which contain A are multiplied 
by the powers 7 + 0 - a and 7 + a - 2(3 of (t s - 1 ) respectively. Since o = | (3 from (2.6) these 

terms are in balance and provide no new information. An expression for 7 is next found by 

making a balance of all four leading order terms in (2.8). This gives 

7 = 2/?-l, (2.9) 
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and the kinematic condition becomes to leading order 

-i //t"=-|/j/ + /3f/' + /X (2-10) 

where primes denote ^-derivatives. 

Another equation, as well as determination of the free parameter /? is found from the 
Bernoulli equation (2.1c). Substitution of the ansatz (2.2), (2.4) and retention of the largest 
possible terms as t — ► t s -, gives the following equation: 

- 7 (<. - ty-'A + {t. - ty-'ftAt + ( ts - ( 2 . 11 ) 

Wc arrive at (2.11) above by dropping any terms that are smaller than the ones retained 
using the fact that cit > ft alone. Next, we balance the different terms in the equation. The 
terms of orders ( t s - f) 7-1 and (t s - t ) 2 ^ 2/? are in balance by virtue of (2.9) above. The 
surface tension term containing / -1 is chosen to balance the other terms and produces a 
coupling with (2.10) between inertial and capillary forces. This is possible if 

7 — 1 = —a. 


Use of (2.6) and (2.9) that provide a and 7 in terms of ft yields a value for ft and hence a 
and 7 follow, 

( 2 . 12 ) 


„_ 6 4 _ I 

® — y j ft ij > 7 75 


while equation (2.11) becomes 


-\a + 4 -(A ( + '-A\ = -J. (2-13) 

Equations (2.10) and (2.13) form a coupled system which must in general be addressed 
numerically. Before presentation of such results we cast (2.10), (2.13) into a form which is 
readily comparable with the one-parameter family equations of PO. We define G = A^ and 
/ = y/F, differentiate (2.13) with respect to £ to yield the system 


-yF + + (GF)' = 0, (2.14a) 

-G + Ug 1 + \(G 2 y = -(E" 1/2 )'. (2.146) 

Equations (2.14a,b) are exactly the equations found from slender jet theory corresponding 
to ft = |. In the absence of surface tension the term on the right hand side of (2.14b) is 
dropped. As shown below, zero surface tension solutions are possible even though the flow is 
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linearly stable. Before presentation of numerical solutions we give the asymptotic behavior 
of (2.14a,b) for large |£| (for more details see PO). In the absence of surface tension we find 

F~£ 3 G~t* £->oo, (2.15a) 

F~(-0 -1 G £ > oo. (2.156) 

Note that the shapes predicted for zero surface tension grow on one side (as £ — *• oo here) 
and decay at the other end of the domain. This is in contrast to the surface tension case as 
the following asymptotic forms indicate, 

F ~ £ 3 G £ — *■ oo, (2.16a) 

F~(-£) 3 G~(-0"* -oo. (2.166) 

Numerical solutions of (2.14a,b) showing the variation of F and G with £ are given in 
Figures 1 and 2 for zero surface tension. Verification of the asymptotic forms (2.15a,b) is 
included in Figures 3 and 4 for the behavior of F and G respectively. This is achieved by 
constructing logarithmic plots of the numerical solutions for large positive and negative £. 
The slopes of the lines so obtained were computed from the data by a least squares fit and are 
indicated on the Figures. It is seen that agreement is very good and improves as the range of 
integration is extended. The results for non-zero surface tension are presented in Figures 5 
and 6 with the corresponding logarithmic plots confirming the asymptotic behavior in Figures 
7 and 8. The main qualitative difference between the two cases is that in the absence of 
surface tension the shape is monotonic tending to zero at one end of the similarity region 
and growing algebraically at the other end, while in the presence of surface tension the shape 
function grows algebraically (at the same asymptotic rate) at both ends of the domain. The 
former type of behavior is also obtained for highly viscous jets analyzed in Section 3. 

2.1 Conservation integrals at breakup. 

The inviscid jet flow is conservative and solutions to the problem are required to conserve 
the total momentum and energy as well as mass. It can be seen from the solutions just 
constructed that velocities are becoming unbounded as the jet radius shrinks to zero. It is 
important, therefore, to verify the consistency of these solutions with conserved quantities of 
the system. We define the total mass, momentum and energy of an axially periodic portion 
of the jet by 7i, I \ and / 3 respectively. These are given by 

,D ,5 

l\ = 2irp j J rdrdz , (2-17) 
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( 2 . 18 ) 


Ij = 2ir pi [ l<t\ + <p‘,)' r ‘rdrdz, 

J-D JO 

/ 3 = 27 Tp f f + 4> 2 z )rdvdz + 27ra f Sdz , (2.18) 

J-D Jo 2 J-D 

where <r is the surface tension coefficient and the quantities (2.17)-(2.19) are in dimensional 
forms. The contributions of the integrals I u I 2 and / 3 in the similarity region are found 
by substitution of the leading order solutions from the ansatz (2.2a-c) to give the following 
quantities 

7i~ (*.-*)“ I" Fd(, C 2 - 19 ) 

J — oo 

h ~ (t, - 1)“ /” FGdt, (2.20) 

h ~ ((, - tf [°° FG'd(. (2.21) 

The contributions to the energy integral (2.18) of the kinetic and surface energy in the 
similarity region are of the same order and form and both appear in (2.21); the expression, 
therefore, is valid both in the absence and presence of surface tension. We have shown 
analytically and numerically that F and G become unbounded as |£| —* oo. In order to 
evaluate the contributions (2.19)-(2.21) on the large scale flow we need to transform to outer 
variables. The main contributions to the integrals come from the asymptotic forms (2.1oa,b) 
for zero surface tension or (2.16a,b) for non-zero tension. The behavior at infinity, however, 
is the same for both cases and the results that follow are equally valid. Consideration of the 
leading order contributions of (2. 19)-(2.21) gives 

A ~ (t* - t)^C 4 = z 4 , ( 2 - 22 ) 

/ 2 ~(t s -t)T^ = (t s -t)i^, (2.23) 

= < 2 - 24 ) 

We can conclude from (2.22)-(2.24) that as the singular time is approached, the solutions 
constructed here are consistent with the conservation integrals of the problem. More specif- 
ically the total mass, momentum and energy in the breakup region tends to zero as the 
singular time is reached. This theoretical result can be of used as a criterion in numeri- 
cal simulations, for example, since near the singular time the total energy of the system is 

concentrated away from the similarity region. 

3 Theory for viscous jets 

In this section we consider the collapse of a viscous thread of fluid of undisturbed radius R 
under capillary instability. In numerous applications characteristic Reynolds numbers are 
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small and the flow is governed by the Stokes equations. Such flows are driven by capillary 
forces which in turn provide a scale for the flow velocities. In non-dimensionalizing the 
equations, therefore, the following scales are used: 


(r,z) = R(r,z) (u,w) = ~{u,w) p = —p t = — t, 

it 7 


where cr is the surface tension coefficient and // the fluid viscosity. The non- 
equations and interfacial boundary conditions become (dropping the bars): 

-dimensional 

A 1 

All ~U = p r , 

V* 

(3.1a) 

Aw = p z , 

(3.16) 

- ( ru ) + w z = 0, 

r 

(3.1c) 

d 2 Id d 2 

dr 2 + r dr + dz 2 

On r = S(t, z ) we have 


( u z + in r )(l — Si) + 2 u r S z — 2 w z S z = 0, 

(3. Id) 

p - 2u t - (-p + 2w z )S 2 z + 2(u z + w r )S z = 


- (s„ - i(l + Si)) (1 + 

(3.1e) 

u = S t + wS z . 

(3.1/) 


The interfacial conditions (3.1 d-f ) represent the tangential stress balance, normal stress 
balance and the kinematic condition respectively. An additional condition is regularity of 
the velocity field on the jet axis r = 0. The problem (3.1a-f) poses a formidable analytical 
task. Analytical studies are usually confined to linear stability and in the fully nonlinear 
regime the problem has been addressed numerically (see Introduction). In order to draw 
up an analogy with the inviscid analysis of Section 2, we can introduce a streamfunction ip 
defined by u = —\ip z , w = }ipr so that the continuity equation (3.1c) is satisfied. Elimination 
of p between (3.1a,b) yields a single equation for ip, namely 

E 4 iP = 0, E 2 = A - -S-. (3.2) 

r or 

The flow field in the viscous case is governed by a fourth order equation as opposed to the 
second order Laplace equation for the potential in inviscid flows. As seen below the analysis 
of viscous jet pinching is slightly more involved. In what follows we choose to work with 
primitive variables rather than with the streamfunction ip. 
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Following the ideas developed in Section 2, we look for singular terminal states of (3.1a-f) 
according to the ansatz 

r = (t. - t)^ z = (t s - S = (t, - 0V(O. 

W = (t 3 -tyw{t,y, 0 u = (t s -ty +a - f) U(t,y,0 p=(t s -t)~ a P. (3.3) 

The expression for u follows from the continuity equation (3.1c) once w is specified, and the 
scaling for the pressure p is a consequence of the normal stress balance equation (3.1e) since 
capillary instability drives the dynamics. As before, the geometry of pinching dictates the 
inequality 

a > j3. (3-4) 

Consideration of equation (3.1b) for w indicates that if a leading order balance is made 
between the radial derivative terms of w and the pressure gradient term (this requires 7 = 
a — the following leading order solution for W arises 

W = + A{Q, (3-5) 

where A(() is an unknown function. Ail inconsistency appears now if (3.5) is substituted 
into the tangential stress balance equation (3. Id). The leading order contribution to (3. Id) 
after substitution of the scalings (3.3) is simply W y = 0 on y = /(£)■ Application of this 
condition to the solution (3.5) found above implies that = 0 leading to an inconsistent 

solution P = const.. A consistent solution is constructed by expanding W (and therefore U) 
in appropriate powers of ( t s — t) with the leading order contribution of W being independent 
of y so that the tangential stress balance is satisfied identically to leading order. Higher order 
corrections from (3. Id) then enter to force a nontrivial solution. The appropriate expansions 


for w and u become, then, \ 

W = (t s - ty (w 0 (0 + (i ts - tfm + (ts - tyw 2 + ...), (3.6a) i 

« = (ts - ty +a - p (Uo(y, 0 + (ts - tfUr + (t, - tyu 2 + ...). (3.66) : 

The solution for U 0 which is regular at y = 0 follows from (3.1c) I 

Uo = -\yW oi . (3.7) - 

Next we obtain the solutions for W\ and U\. Using (3.3) and (3.6a) in equation (3.1b) gives 

<«. - <r- 2 " (|; + i J) ((*. - o-tv, + «. - tyw, + . . .) ; 

+(t s - ty-^Woa + ... = (t s - t)- a - 0 Pt. (3.8) 
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Balance of the leading order terms involving Wo, Wj and P gives 

a = \ & i = . 


(3.9) 


and a solution for W\ follows by integration with respect to y since the pressure is independent 
of y by the normal stress balance (see comments below also) 

"'i = - W m ) + A(f), (3.10) 

where A is some function of £. Using (3.10) and (3.1c) gives 

Vi = -—»»(/>« - Woe«) + \yA t . (3.11) 

Implicit in the solutions (3.10), (3.11) is the assumption that P y = 0. Using the values (3.9) 
for a and 7 shows that the leading order balance in (3.1a) is 


' \_d_ 
,dy 2 + ydy 



(3.12) 


Now substitution of the solution (3.7) for Uq into (3.12) shows that the left hand side is zero 
and so 

Py = 0 . (3.13) 

Equation (3.13) is crucial for the consistency of the theory. 

Next we consider the tangential stress balance. Using the values (3.9) gives the following 
equation to leading order 


Uoi + w ly + 2Uoyf - 2 Wotf = 0 y = f((). (3.14) 

Use of solutions (3.7) and (3.10) into (3.14) gives an ordinary differential equation for Wo(£), 

= (3.15) 

An expression for P in terms of Wo is available from the normal stress balance equation. It 
is found that to leading order several terms in this equation are in balance (according to the 
previously found values (3.9)) the result being 

P- 2!/„» = j y = /(£). 

Substitution of this into (3.15) and integration gives the solution 

=— + •£, (3.16) 
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where k is a constant of integration. With Wo known in terms of /, an equation for / arises 
from the kinematic condition (3. If) and use of the ansatz (3.3) along with the appropriate 
derivative transformations (see Section 2 also). To leading order, therefore, the kinematic 
condition becomes 

(t s - ty +a - 0 Uo = (ts - + P£f) + (t, - tr+^Wof. 

A balance of terms is achieved if 7 + a — /3 = a — 1 which along with (3.9) gives 

a=l /J = | i = -\. (3.17) 

The kinematic equation is 

(ni + ft)/+(- I + 5"S)/-o, 

which on elimination of Wq from (3.16) becomes 

( w> + h) f + ( _/ _ g + h ) = °' (3 - 18) 

Equations (3.16) and (3.18), therefore, provide a coupled system of first order equations to 
be solved subject to conditions at £ = 0, for instance. In subsections 3.1 and 3.2 that follow 
we consider analytical and computational aspects of this system. 

First, some comments on the consistency of our theory with the original assumption that 
the flow is in the Stokes regime are necessary. It is seen from the values (3.17) that the 
axial and radial velocities (see (3.6a,b)) have size 0((f s — t) -1 ^ 3 ) and 0(1) respectively at 
the time of pinching. The magnitude of the velocity vector is becoming infinite, therefore, 
and the possibility of violation of the Stokes flow assumption arises. The terms ignored are 
the unsteady and convective terms of the Navier-Stokes equations which read 

R e (u t + uu r + uw z ) = — p r + An — (NS 1) 

R e (w t + uw T + ww 2 ) = —p z T Aw, (N 5'2) 

where R e = is the Reynolds number. In deriving the similarity equations (3.16), (3.18) we 
use information from IJq,U\ in (NS1) and W$,W\ in (NS2). The smallest terms considered, 
therefore, are of 0((t s - t)~ 4/3 ) in (NS1) and 0((t s - t)~ 5 ^) in (NS2), coming from the 
Laplacian of U\ and Wi respectively. The largest possible terms which would enter through 
the unsteady and convective terms on the left hand side of (NS1) and (NS2) are easily 
calculated to be of orders 0((t s - <) -1 ) and 0((t. s - t)~ 4/3 ) respectively. As t — » t. s ~, 
therefore, the validity of the Stokes approximation and the ansatz leading to (3.16), (3.18) 
is consistent; higher order terms eventually enter via unsteadiness and nonlinearity in (NS1) 
and (NS2) but are of no consequence to the leading order balances constructed here. 
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3.1 Analysis of the similarity equations 

Elimination of Wq between (3.16) and (3.18) yields a single second order equation for / 


which can be written as: 


r+i 

6 2) it* 


P + f- 


(3-19) 


At first sight, equation (3.19) appears to admit a class of even solutions which, without loss 
of generality, may satisfy /( 0) = /o, f'( 0) = 0. A series solution for small however, with 
the above initial conditions is found to yield a constant solution /(£) = / 0 . For general 
values of £ it can be shown that acceptable solutions to (3.19) produce monotonic functions. 
Two distinct cases need to be considered, (a) k < 0, (b)k > 0. 

We consider the case k < 0 first. If we assume that /' vanishes at some finite value 
of £ = 6) say, then by (3.18) we deduce that W^(£ 0 ) = 2. Using this in (3.16), however, 
requires = 2 which is impossible since k <0 and /(£ 0 ) > 0, showing that f 

is non- vanishing and so /(£) is monotonic. This result can also be proved by construction 
of closed form expressions for /'. We integrate (3.19) once with respect to £ so that 

| = AGU), 

if™ 0 ™ - Wf H0ry ^ 

where A is a constant. Equation (3.20) can be integrated explicitly. There are three distinct 
cases: (i) k — 0, (ii) k < — (iii) — ~ < k < 0, giving the following solutions. 

Case (i), k = 0 

r*t t\ _ J n, \ 


G(f) = 


(/ + 1 / 6 ) 


(3.21a) 


Case (ii), k < -1/72 (k = -2A 2 ) 


G{f) = 


(P + |/ + A 2 )3 


V A2 - 15 


f+A 


(3.216) 


Case (iii), -1/72 < k < 0 (jt = (t + 




-1+ J- 

3 ~ 9 /* 


(3.21c) 


These expressions show that G(f), and hence /', is non-vanishing; cases (i) and (ii) are 
straightforward while in case (iii) the observation |//| < -jb yields the desired conclusion. 
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When k > 0 the expression for G is identical to (3. 21c) with fi as given there. It follows 
that for k > 0, 

h <ll< °°' 

with equality on the lower bound already covered by the case k = 0 in (3.21a). Using this ^ 

inequality we see from (3.21c.) that at positions where f = fi — the derivative / can (i) 
become infinite, (ii) be finite and non-zero, (iii) become zero, depending on the value of the 
exponent — ^ + <j~. In terms of /i these three cases are 

1 1 * 

">3 ’ " = 3 ’ l ‘ < r 

respectively. In what follows we show (by sketching the behavior rather than presenting 
rigorous proofs) that case (ii) alone allows admissible solutions. 

We consider the range | < /z < oo first. As mentioned above the slope /' becomes 
infinite at positions £ 0 given implicitly by /(£ 0 ) = V ~ n- The i nter facial slope becomes 
double valued; a local analysis in the vicinity of /o = ft — gi ves 

We note that an infinite slope in / violates the ansatz established by the scalings (3.17) since 
spatial variations in the axial direction become much larger than the variation in the radial 
direction. 

Next, we consider the values — < n < It is seen from (3.21c) that the derivative 
vanishes when f = \2 — M since the exponent of the last factor satisfies 0 < 3 9^ ^ 

This inequality for the exponent, however, implies that higher derivatives of the interfacial 
amplitude in the similarity region are unbounded and this case is dropped since we are 
interested in smooth solutions. 

The final possibility is // = | (k = ^) which casts (3.21c) into 

(?(/) = /(/ + ( 3 - 21<i > 

Equations (3.21a) and (3.21d) can be integrated to yield closed form expressions for £ in 
terms of /; the solutions for / are implicit, therefore, and are not very illuminating besides 
the checks they provide for the numerical and asymptotic work. These solutions are 
fe = 0, (3.21a) 

+ §/(/; -ij) = + B°, (3.22a) 

h=h (3.21 d) 
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ill (Ml I I IP MW! If III ll ill iWriMIMMIIPIIlH 


where 


(3.226) 


5 l/3 

Hh wr3 ) = M + B u 


-K 2 ln 


/(/; *) = \u + « 3 )* + K*ln ((/ + « 3 )i - «) 

(( /+K 3)| + K(/ + y£ 3) i + Ki ) + W taii _ ] 


The constants /t 0 , A] , -Bo, Si are found from the initial conditions on / and /' at £ — 0 for 
instance. The behavior of the solutions (3.22a, b) for large |£| are easily obtainable. Without 
loss of generality we assume A 0 ,i > 0 to obtain 

k = 0 

'2y4 0 \ 3/i 


/~(^)> <-*00, 

, 1 1 , 
f 7T7 

03 4 


(3.23a, 6) 


fc — I 

K ~ 3 


f 




Z -> oo, 


y4i 12 1 / 3 

/ ~ exp(- 5l/3 0 £ -4 -oo. (3.24a, 6) 

For general values of k the asymptotic behavior is more easily obtainable directly from 
equation (3.19). Assuming that f'( 0) > 0, we look for solutions that grow as £ — > -f-oo but 
decay as £ -» -oo (such behavior is supported by the closed form solutions (3.22a, b) as well 
as numerical calculations). To leading order, equation (3.19) takes the forms 


f" 1 /' , 

7 = 5J />>1 ' 

f r r" f 

J = J (*/») . J = 2 J (* = «) /« 1 . 


(3.25a) 


(3.256, c) 


It can be seen that for large positive ( the asymptotic form is independent of k , in agreement 
with (3.23a) (3.24a) above. The asymptotic forms are, therefore, 

/ ~ f /2 oo, (3.26a) 

/ ~ exp(- 9 |£|) (k^O) / ~ l^r 1 (A: = 0) ( -> -oo, (3.266) 

where q is a positive constant which can be determined numerically. Given the behavior of 
/ the asymptotic forms for the axial jet velocity Wo follow from equation (3.16), 

Wo ~ const. + 0(£ _1/2 ) £ -4 oo, (3.27a) 
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Wo ~ exp(2g|{|) (k ± 0) W 0 ~|£| 3 (k = 0) -oo. (3.276) 

In the following section we present numerical solutions of the similarity equations with 
particular emphasis on the verification of the analytical findings of this section. 

3.2 Numerical solutions of the similarity equations 

Equations (3.16), (3.18) or equivalently equation (3.19) were solved numerically by specifying 
conditions at ( = 0 and integrating out to -f oo and -oo. Typically we worked with the second 
order equation (3.19) with initial conditions /( 0) = /o and /'( 0) = /i > 0. Integration to an 
appropriately large positive and negative value of £ was achieved by a fourth order Runge- 
Kutta method. Representative solutions are considered for the three cases (l) k < 0, (u) 
k = 0, (iii) k = §. We note that the distinctions between different negative values of k are 

not important in the numerical calculations. 

Figure 9 shows the results for / and W 0 of a case having k = -1 and initial conditions 
f(0) = 1 and /'(0) = 1. To establish the asymptotic behavior of the solutions for large 
|£| we construct logarithmic plots of the relevant parts of the solutions. Results for large 
negative £ are shown in Figure 10 which shows the variation of ln(f) and ln(w) with £. The 
computed slopes included on the curves were calculated by a least squares fitting of the data. 
According to the asymptotic forms (3.26b) and (3.27b), w grows exponentially at twice the 
rate at which / decays exponentially. The computed slopes of 2.269 and —1.138 for hi(w) 
and ln(f) respectively are in full agreement with this finding. The asymptotic form (3.26a) 
for large positive £ is confirmed in Figure 11 with very good agreement. Verification of the 
behavior of W 0 was also obtained by plotting the variation of W 0 with £ -1/2 for large 

For k = 0 and the same initial conditions as above qualitatively similar solutions are 
obtained for both / and W 0 and representative solutions are given in Figure 12. A comparison 
between the numerical and asymptotic solutions for large positive £ according to (3.23a) is 
described next. The behavior (3.23a) is equivalent to 

<»(/) ~ l‘n (j) + j/n (|) + f '"«)• 

The value of the intercept follows from (3.20) when the values /( 0) = /'( 0) = 1 are used, 
and is equal to -0.224. A logarithmic plot of the numerical solution produced a slope of 
1.498 and an intercept -0.208 (these results were obtained by integration to £ = 500; these 
values improve as the maximum integration range is increased). The results for k = 3 are 
qualitatively similar to those described above and are not included here. 
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4 Conclusions 


A formal theory has been developed to describe the breakup process of highly viscous or 
inviscid jets in air. The similarity solutions obtained here constitute inner solutions in the 
vicinity of regions where the jet radius is going to zero which when matched with outer solu- 
tions of the Stokes or Euler equations respectively, provide a full description of the breakup 
phenomenon. In general the outer solutions are fully nonlinear and must be computed numer- 
ically. In numerical calculations of jet breakup (by boundary integral methods for example), 
it is important to have a rational way of continuing the computations beyond the change in 
topology necessitated by the physics. The similarity solutions presented here provide such 
a possibility since accurate initial conditions for the pinched part of the evolution, consis- 
tent with the full equations of motion can be easily calculated. Direct simulations near the 
pinching times contain stiffness problems which are likely to produce inaccuracies. Another 
example of the need of accurate initial conditions for the jet evolution beyond breakup is the 
experimentally and computationally observed behavior of a residual fluid tube at breakup 
being accelerated by surface tension to form a jet which can go through the mother drop 
it was connected to. The speeds of such jets depend on interfacial tension forces which are 
proportional to the interfacial surface area and hence the shape at breakup. 

As an illustration of the mechanics of the matching process we consider the results of 
Section 3 for highly viscous jets. Matching is achieved as |£| —► oo and so the asymptotic 
forms of the solutions are essential in achieving this. Using (3.26a), for example, we have 
/ ~ £ 3 ^ 2 as £ * oo. Changing to physical (outer) variables using the ansatz (3.3) and the 

scales (3.17) gives the following condition to be satisfied by the outer solution as z — * 0 (we 
have assumed, without loss of generality that pinching occurs at z = 0), 

. S~z*. 


The multiplicative constant that goes along with this behavior and which can be found 
from the outer solution, fixes the inner problem whose solution is easily obtained from the 
similarity equations. 

We have also found a qualitative difference between the breakup forms of inviscid jets 
with surface tension and their highly viscous counterparts. In the former case our solutions 
indicate that at breakup the interface is more or less symmetric with growth to mother drops 
on either side of the pinch point. In the latter flows (as well as inviscid jets with zero surface 
tension), our similarity solutions indicate that at breakup there are thin, and relatively long, 
fluid tubes which taper to zero diameter on one side and grow to connect to a mother drop 
on the other side. These findings are in full qualitative agreement with photographs from 


17 



various experiments on both highly viscous jets (Tjahjadi et al. 1992) as well as inviscid 
ones (Donnelly and C.laberson 1966, Chaudhary and Maxworthy 1980a, b). 
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Zero surface tension (inviscid) 
f(0)-1 g(0)-l beta*4/7 



5 

• Figure 1 Inviscid jet with zero surface tension; similarity solution /(£). 


Zarosurtaca tension (inviscid) 
f(0)-1 o(0>— i bata-4/7 



• Figure 2 Inviscid jet with zero surface tension; similarity solution </(£). 
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• Figure 9 Highly viscous jet with surface tension. Similarity solutions W 0 (£) /(£)• 
k = — 1, initial conditions /( 0) = 1, /'( 0) = 1. 


Behavior at minus infinity 
k— 1 1(0)- 1 niH 



• Figure 10 Highly viscous jet with surface tension. Asymptotic behavior as £ — * +oo. 
k = — 1, initial conditions /(0) = 1, /'( 0) = 1. 
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Figure 12 Highly viscous jet with surface tension. Similarity solutions Wo(£) /(£) 
k = 0, initial conditions /( 0) = 1, /'( 0) = 1. 
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